function hfile = srf_CurvatureColor(hfile, csmooth)
% SRF::CurvatureColor  - color according to curvature
%
% FORMAT:       srf.CurvatuteColor([csmooth])
%
% Input fields:
%
%       csmooth     amount of smoothing (nr of neighbors to mean)
%
% No output fields.

% Version:  v0.7b
% Build:    7090215
% Date:     Sep-02 2007, 3:51 PM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'srf')
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if nargin < 2 || ...
   ~isnumeric(csmooth) || ...
    isempty(csmooth) || ...
    isinf(csmooth(1)) || ...
    csmooth(1) > 5
    csmooth = 0;
end

% build curvature
p = bc.VertexCoordinate;
n = bc.VertexNormal;
nei = bc.Neighbors;
numv = bc.NrOfVertices;
curve = zeros(numv, 1);

% iterate over vertices
for vc = 1:numv
    
    % set initial curvature estimate to 0
    cp = 0;
    
    % get reference point
    rp = p(vc, :);
    
    % get neighbors
    pcnx = nei{vc, 2};
    
    % iterate over neighbors
    for pc = pcnx(:)'
        
        % get vector and normal
        tv = p(pc, :) - rp;
        nv = n(pc, :);
        
        % add to curvature estimate
        cp = cp + acos((tv * nv') / sqrt((tv * tv') * (nv * nv')));
    end
    
    % set into array
    curve(vc) = cp;
end

% apply some smoothing
for cc = 1:csmooth
    
    % reduce by mean
    rcurve = curve - mean(curve);
    
    % fill array anew
    curve = zeros(numv, 1);
    for vc = 1:numv
        
        % get curvature estimate
        cp = rcurve(vc);
        
        % get neighbors
        pcnx = nei{vc, 2};
        
        % add curvature of neighbors
        for pc1 = pcnx(:)'
            cp = cp + rcurve(vc);
        end
        
        % and store back
        curve(vc) = cp;
    end
    
    % limit data
    cstd  = std(curve);
    cmean = mean(curve);
    curve(curve < (cmean - cstd)) = cmean - cstd;
    curve(curve > (cmean + cstd)) = cmean + cstd;
    curve = fix(0.5 + (curve - min(curve)) / (max(curve) - min(curve) + eps));
end

% split into 0 and 1
mc = mean(curve);
curve = double(curve(:) >= mc);

% create color indexing
curve(:, 2:4) = 0;

% set into file
bc.VertexColor = curve;
bvqxfile_setcont(hfile.L, bc);
